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Abstract 

A bivariate ensemble model output statistics (EMOS) technique for the postprocessing of 
ensemble forecasts of two-dimensional wind vectors is proposed, where the postprocessed 
probabilistic forecast takes the form of a bivariate normal probability density function. The 
postprocessed means and variances of the wind vector components are linearly bias-corrected 
versions of the ensemble means and ensemble variances, respectively, and the conditional cor- 
relation between the wind components is represented by a trigonometric function of the en- 
semble mean wind direction. In a case study on 48-hour forecasts of wind vectors over the 
North American Pacific Northwest with the University of Washington Mesoscale Ensemble, 
the bivariate EMOS density forecasts were calibrated and sharp, and showed considerable im- 
provement over the raw ensemble and reference forecasts, including ensemble copula coupling. 

1 Introduction 

The past two decades have seen a change of paradigms in weather forec asting, in that ensemble 
predi ction systems have been developed and implemented operationally (|Leutbecher and Palmer , 
|2008|) . Ensemble systems seek to reflect and quantify sources of uncertainty in numerical weather 
forecasts, such as imperfections in initial conditions and incomplete mathem atical representations 
of the atmosphere. Despite th e ubiquitous positive spread-skill relationship (|Whitaker and Loughe . 



1998; iGrimit and Mas si. 120021). ensembl e forecasts tend to be biased, and typically they are un- 



derdispersed (|Hamill and Coluccil |1997|) . in that the ensemble spread is too small to be realistic. 
Furthermore, differing spatial resolutions of the forecast grid and the observation network may 
need to be reconciled. 

To address these shortcomings, various techniques for th e statistical postprocessing of ensemble 
model output have been developed (|Wilks and HamillL 120071). with ensemble model output statistics 
(EM OS) or nonhomogeneous Gaussian regression (|Gneiting et aUl2005l ; lThorarinsdottir and GneitingL 
|2010|) being a state of the art method. The EMOS technique transforms a raw ensemble forecast 
into a predictive probability density function, and simultaneously corrects for biases and dispersion 
errors . EMOS method s have been developed for temperature and surface pressure (|Gneiting et al. , 



2005; lHagedorn et all 120081 : iKann et all 120091) . where the predictive density is normal and the 



method is oft en referred to as non homogeneous Gaussian regression, for quantitative precipitation 
(Wilksl 120091) . and for wind speed (IThorarinsdottir and Gneitingll2010l ; lThorarinsdottir and Johnsonl 
201 1|) . In all these implementations, the predictive density applies to a univariate weather quantity. 
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Figure 1: Raw ensemble and EMOS postprocessed forecasts of surface wind vectors at stations in the 
Olympic Peninsula and Puget Sound area in the US state of Washington, valid October 20, 2008 at 00 
UTC, at a prediction horizon of 48 hours. The eight members of the University of Washington Mesoscale 
Ensemble (UWME; Eckel and Mass 2005) are shown as gray dots. The 75% prediction ellipse and the mean 
vector for the EMOS density forecast are shown in dark gray, and the verifying wind vector is represented 
by a small black circle. 



In this current paper, we propose and develop an EMOS technique for a bivariate weather 
quantity, namely surface wind vectors, comprising both zonal and meridional components or, in 
an alternative but mathematically equivalent representation, wind speed and wind direction. Prob- 
abilistic forecasts of wind conditions are critical in a wide range of applications, including air 
traffic control, ship routing, recrea tional and competitiv e sailing, and wind energy, where their so- 
cietal and monetary value is huge ([Marquis et all 1201 lb . Until very recently, wind speed and wind 
direction have been a ddressed independently in statist i cal postprocessing, without taking depen- 
dencies into account (IBao et all 12010c ISloughter et all 120101 : iThorarinsdottir and Gneiting , l20ld : 



Thorarinsdottir and JohnsonLl201 lh . However, in many of the aforementioned applications it is im- 



portant to honor the full information about the bivariate structure of the future wind vector that is 
provided by the ensemble. Thus, our EMOS postprocessed forecasts take the form of elliptically 
symmetric bivariate normal densities, as illustrated in Figure Q] in an application to the University 
of Washington Mesoscale Ensemble (UWME; Eckel and Mass 2005). 

The remainder of the paper is organized as follows. In Section |2] we provide the details of 
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the bivariate EMOS technique. A case study is presented in Section [3} where we consider 48- 
hour ahead forecasts of wind vectors over the North American Pacific Northwest in 2008 based on 
the eight-member UWME. The paper closes in Section @] where we hint at future developments 
and discuss the s i milarities and di f ferenc es between our EMOS technique, the BMA approach of 
Sloughteil d2009h : Isioughter et~aD d201lb. ensemble co pula coupling (ECC; Schefzik 2011) and 



the postprocessing method proposed by iPinsonl (|201ll) . all of which are directed at the bivariate 
postprocessing of ensemble forecasts of wind vectors. The Appendix describes our verification 
methods. 



2 Ensemble model output statistics for wind vectors 

A wind vector is determined by wind speed and wind direction, or by its zonal (west-east) and 
meridional (north-south) components, which we denote by u and v, respectively. We now develop 
an ensemble model output statistics (EMOS) method for wind vectors, where the postprocessed 
probabilistic forecast takes the form of a bivariate probability density function. The method is 
tailored to ensembles with relatively few members, such as the eight-member University of Wash- 
ington Mesoscale Ensemble (UWME; Eckel and Mass 2005), and we illustrate it using forecast 
and observation data from this ensemble. 



2.1 Bivariate normal distribution 



Our EMOS postprocessed forecast takes the form of an elliptically symmetric, bivariate normal 
probability density function for the wind vector (u, v), with the parameters of this distribution being 
specified in terms of the ensemble forecast. The analytic form of the bivariate normal probability 
density function is 
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The task now is to specify the five parameters in CD, namely the mean values, p, u and p, v , of the wind 
vector components u and v, the corresponding marginal variances, a 2 and a 2 , respectively, and the 
correlation coefficient, p uv , between the wind components, in their dependence on the ensemble 
forecast. 



Bivariate normal density forecasts for wind vectors have al so been propos ed by iGneiting et al 



(|2008|) . though in very crude form. The ingenious method of IPinsonl (|201ll) estimates a dilation 
and translation of an ensemble forecast of wind vectors based on bivariate normal densities, and 
our approach is very similar in its tre atment of the m ean and variance parameters. However, major 
differences between the approach of IPinsonl (|201l|) and our method lie in the form of the post- 
processed forecast, which is a probability density function in our case, rather than a dilated and 
translated ensemble, and in the explicit modeling of the correlation coefficient p uv in our method. 
For a more detailed comparison and recommendations for a judicious choice of the most appropri- 
ate method, given any particular ensemble and task at hand, we refer to Section HI 
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In the description that follows, we consider a general ensemble with m members, and denote 
the individual wind vector forecasts by (ui, v{), . . . , (u m , v m ), respectively. 



2.2 Means 

In our standard implementation of the bivariate EMOS technique, the means \i u and p v are bias- 
corrected versions of the respective ensemble means, in that 

fJ'u — a u + b u u and //„ = a v + b v v, (2) 

where u = — Y^iLi u i an d v — — YlT=i v i- The bias correction parameters a u , b u , a v and b v are es- 
timated from training data, in ways described below. In a slightly more ambitious implementation, 
the mean components fi u and fj, v are affine functions of the individual ensemble member forecasts, 
namely 

Vu = a u + b u ,iui H h b u>m u m and [i v = a v + b Vf iVi H h b Vym v mi 

where the regression parameters a u , b Ut i, . . . , b u>m , a v and b v ,\ , . . . , b v>m are estimated from training 
data. This general version applies to ensembles with non-exchangeable members only and reduces 
to the standard version when b Ut i = ■ ■ ■ = b Ujm = — and b V) i = ■ ■ ■ = b v>m = — . In our experiences 
with the UWME, which has non-exchangeable members, the general version gave only very slightly 
improved predictive performance, and so we report results for the standard implementation © only. 



2.3 Variances 

We specify the marginal variances a\ and a 2 of the bivariate normal density forecast (OQ) as affine 
functions of the respective ensemble variances, in that 

a l = °u + d u sl and a 2 v = c v + d v s 2 v: (3) 

where s 2 u = — Y^T=i( u i ~ an( ^ s 1 = ~ Y^T=i( v i ~ ^) 2 - ^he dispersion correction parameters c u , 
d u , c v and d v are estimated from training data, as described below. To guarantee the nonnegativity 
of the variances, we constrain the parameters to be nonnegative, using the technique described by 



Thorarinsdottir and Gneitingl (|2010|) 



2.4 Correlation coefficient 

The key characteristic and major innovation in our work is the explicit modeling of the correlation 
coefficient p uv of the postprocessed bivariate normal density forecast CD- 

To motivate our specification of p uv , we consider ensemble forecast and observation data from 
the UWME in calendar year 2007, comprising a total of 23,250 forecasts cases at 79 meteorolog- 
ical stations in the Pacific Northwest. Figure [2] distinguishes nine sectors for (u,v) wind vector 
forecasts and observations. Figure [3] shows wind vector observations in 2007, conditionally on the 
ensemble mean wind vector falling into any given sector. The conditional distributions are partly 
elliptically contoured, particularly in the first sector, and partly skewed, with an orientation that 
depends strongly on the sector, and they reflect the discretized nature of wind observations, as 
discussed in more detail in Section [31 
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Figure 2: Sectors for wind vector forecasts and observations in the (u, v) plane. Sector 1 is the circular 
region that is centered at the origin and corresponds to wind speeds less than or equal to two meters per 
second. Sectors 2-9 are assigned clockwise, with Sectors 2-3 corresponding to a south-westerly, Sectors 
4-5 a north-westerly, Sectors 6-7 a north-easterly, and Sectors 8-9 a south-easterly wind. 



The left-hand panel in Figure |4] plots the correlation coefficient between the wind components 
in the scatter plots in Figure [3] as a function of the wind directions that correspond to the centers 
of sectors 2-9. The panel demonstrates that the ensemble mean wind direction ought to have a 
profound influence on the correlation coefficient p uv in the postprocessed EMOS density forecast 
(0Q). Thus, we model the correlation coefficient p uv as a trigonometric function of the ensemble 
mean wind direction, 9, measured in degrees, in that 

p uv = rcos(J^(k6 + (p)j +s, (4) 

where the parameters r, s, k and ip are estimated from training data, in ways described below. The 
coefficients r and s concern the overall magnitude of the correlation coefficient and need to satisfy 
l r l + \s\ < 1. The parameter k corresponds to the number of periods of the trigonometric function, 
which we constrain to be either k = 1, 2 or 3, and the direction ip encodes phase information. 

We apply the correlation model © to forecasts with ensemble mean vectors in all nine wind 
sectors, including the first sector. Alternatively, if the ensemble mean wind vector falls into the 
central first sector, one can take p uv to be equal to the empirical correlation coefficient in the cor- 
responding scatterplot. In our experiences with the UWME, the two approaches resulted in nearly 
identical predictive performance. 



2.5 Estimation 

Our EMOS postprocessed density forecast for a wind vector takes the form of the bivariate normal 
probability density (Q}, where the relationships of p u ,p v ,a\, a 2 v and p uv to the raw ensemble fore- 
cast are specified in equations ©, © an d ©. It remains to estimate the parameters that govern 
these equations, and we address this task in three phases. 

Before describing the phases o f our estimation scheme, it is worth noting that we follow 



Thorarinsdottir and Gneitingl (|2010|) and consider two distinct variants, which we call the Regional 



EMOS and the Local EMOS method, respectively. In the Regional EMOS method, only one set of 
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Figure 3: Wind vector observations over the Pacific Northwest in 2007, conditional on the ensemble mean 
forecast falling into one of the sectors defined in Figure |2l The unit for the wind components is meters per 
second. 
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Figure 4: Left: Correlation coefficient between the wind components in the scatter plots in Figure [3] as a 
function of the wind direction that corresponds to the center of the sector. The correlation coefficients and 
observation counts correspond to sectors 2-9, respectively. The dashed curve shows the correlation model 
© as fitted by the weighted least squares method. Right: Same as left, but considering data from Sea-Tac 
Airport only. 



parameters is estimated and used to produce forecasts over the entire ensemble domain, such as the 
Pacific Northwest for the UWME. Training sets thus comprise data from all stations. In contrast, 
the Local EMOS method uses training data from the station at hand only, and thus obtains a distinct 
set of parameters for each station. 

We now turn to the first phase of our estimation scheme, in which we fit the correlation model 
©. We do this offline, once and for all, based on historic, out-of-sample forecast and observation 
data. Specifically, for the UWME, we use data from calendar year 2007 to form the conditional 
scatterplots in Figure [3] and compute and plot the corresponding empirical correlation coefficients, 
as illustrated in FigureUl We then decide about a suitable value for the number of cycles, k, and fit 
the remaining parameters, r, s and cp, of the correlation model © by a weighted non-linear least 
squares technique, using the R function NLS with the w eights being proportional to the number of 
observations in the sectors (|R Development Core TeamL I201 lb . The use of a weighted least squares 
technique is critical, particularly for the Local EMOS technique, as local wind patterns may result 
in very few observations being available in any given sector. This first phase of the estimation is 
done once and for all, using historic data from 2007, and the fitted correlation model is applied 
throughout calendar year 2008, which we took as our test period. 

The left-hand panel in Figure 0] shows the fitted correlation model for the Regional EMOS 
method, where we chose k = 2 and obtained weighted least squares estimates of r = 0.20, 
s = —0.15 and <p = —61.9 degrees, respectively. As noted, the Regional EMOS method uses 
these parameters throughout the UWME domain and throughout the test period in calendar year 
2008. The right-hand panel shows the Local EMOS correlation model at Seattle-Tacoma (Sea-Tac) 
Airport, where we took k — 1 and obtained weighted least squares estimates of r = 0.24, s = 0.07 
and ip = 70.5 degrees, respectively. The fitted Local EMOS correlation mo dels at the rema ining 
stations in the UWME domain are provided and illustrated in the appendix of lSchuhenl (|201l|) . 

In contrast to the first phase of our estimation scheme, the second and third stages proceed on- 



7 



line, that is, they use rolling training periods consisting of data from the recent past. The Regional 
EMOS method uses all available data from the Pacific Northwest from the last n days prior to the 
forecast being made. Missing data are simply omitted from the training set. For the Local EMOS 
method, the training period comprises data from the station at hand from the n most recent days 
where forecasts and observations were available. In either case, we talk of a sliding n-day training 
period. 

In the second phase of our estimation scheme, the parameters a u , b u , a v and b v in the specifi- 
cation © of the mean vector are estimated from the training data by standard linear least squares 
regression. In the third phase, the parameters c u , d u , c v and d v in the specification © of the 
marginal variances are estimated on the same rolling training period by the maximum likelihood 
technique, with all other parameters being held fixed. In other words, we maximize the logarithm 
of the likelihood function, namely 

l{cui dui Cvj dy) — ^^( x j\ log f ' ^ (c u , d u , c v , d v ) , (5) 

as a function of the parameters c u , d u , c v and d v in the variance model ©, which are all constrained 
to be nonnegative. The sum in the log likelihood function extends over all locations x and times t 
for which there are data in the training set. Any single term of the form /fa**) (c u , d u , c v , d v ) refers to 
the bivariate normal density (Q]) evaluated at the verifying values u = u^'*' and v = i>( x >*), with fi u , 
\i v and p uv set at the numerical values implied by the mean model © and the correlation model ©, 
based on the parameter estimates from the first two phases of our estimation scheme, and putting 
a u = Cu + duSu and al = c v + d v sl^ x,t \ where s^ x '^ and sl^ x,t ' denote the ensemble variances 
at location x and valid time t, respectively. The optimization is performed numerically with the 
OPTIM function in R, using the Broyden-Fletcher-Goldfarb-Shanno algorithm with initial values 
provided by the previous day's estimates. 

The interpretation of the right-hand side of © as a log likelihood function is valid only if the 
forecast errors are independent between times and locations. While this is usually not the case, an 
alternative interpretation as a mean logarithmic score permits us to view the estimates as optimu m 



score estimates, which are tailored to the estimation of forecasting models (|Gneiting et all 12005) 



2.6 Example 

Figure[5]and Tabled] illustrate the postprocessed Local and Regional EMOS density forecasts of the 
surface wind vector at Sea-Tac Airport, valid October 20, 2008 at 00 UTC, at a prediction horizon 
of 48 hours. Thus, the forecast concerns the same valid time and the same prediction horizon as 
the Local EMOS forecasts illustrated in Figure [Q where the station at Sea-Tac Airport is located 
in the south-east Puget Sound area. The postprocessed density forecasts correct for the biases 
and underdispersion in the raw UWME. The Local EMOS forecast retains the positive correlation 
structure in the raw ensemble and is sharper than the Regional EMOS forecast. 



3 Case study: Forecasting surface wind vectors over the Pacific 
Northwest 

We now consider the out-of-sample predictive performance of our two-dimensional EMOS tech- 
nique in a case study for wind vector forecasts over the North American Pacific Northwest in 2008, 
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Figure 5 : Contour plot of the postprocessed Local EMOS (left) and Regional EMOS (right) density forecasts, 
along with the raw ensemble forecast, of the surface wind vector at Sea-Tac Airport, valid October 20, 2008 at 
00 UTC, at a prediction horizon of 48 hours. The eight members of the University of Washington Mesoscale 
Ensemble (UWME; Eckel and Mass 2005) are shown as gray dots, along with the 90% prediction ellipse, 
which is based on a bivariate Gaussian fit to the ensemble values. The 25%, 50%, 75% and 90% prediction 
ellipses and the mean vector for the postprocessed EMOS density forecast are shown in dark gray. The 
verifying wind vector is represented by the small black circle at (u,v) = (1.45, —0.53). The units are in 
meters per second. 
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Table 1 : Predictive means, variances and correlation for the postprocessed Local and Regional EMOS density 
forecasts in Figure [5] of the surface wind vector at Sea-Tac Airport, valid October 20, 2008 at 00 UTC, at a 
prediction horizon of 48 hours. The parameters refer to the general bivariate normal probability density (OQ), 
with the wind components represented in meters per second. 
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based on the University of Washington Mesoscale Ensemble (UWME; Eckel and Mass 2005). We 
compare to the raw ensemble forecast and various reference techniques, such as ensemble copula 
coupling, and assess the performance of the bivariate EMOS technique when forecasts of wind 
speed are desired only. 



3.1 University of Washington Mesoscale Ensemble 

In our case study, the test set consists of forecasts of surface (10 meter) wind vectors based on the 
University of Washington Mesoscale Ensemble (UWME; Eckel and Mass 2005) with valid date in 
calendar year 2008, at a prediction horizon of 48 hours. The UWME is an eight-member multi- 
analysis ensemble then based on the Fifth- Generation Penn State/NCAR Mesoscale Model (MM5) 
with initial and lateral boundary conditions obtained from operational centers around the world. 
The forecasts are made on a 12 km grid and the region covered is the Pacific Northwest region 
of Western North America, including the US states of Washington, Oregon and Idaho, as well as 
the southern part of the Canadian province of British Columbia. The forecasts were bilinearly 
interpolated from the four surrounding grid points to the observation locations and rotated to match 
the true direction at each station. 

Surface wind vector observations were pro vided by the weather observati on stations in the 



Automated Surface Observing System network ([National Weather Servicel.ll 99 8[) . The vector wind 



quantity studied here is horizontal instantaneous surface wind, where 'instantaneous' means that 
the wind was measured and averaged over the last two minutes before the valid time at 00 UTC. 
The wind vector observations were recorded as wind speed and wind direction, where wind speed 
was rounded to the nearest whole knot, where a knot equals 0.5144 meters per second, while 
values below two knots were recorded as zero. The observations ar e thus discret ized, as is easily 



recognizable in Figure [3] Quality control procedures as described in LBaarsI (|2005|) were applied to 
the entire data set, removing dates and locations with any missing forecasts or observations. 

For calendar year 2008, 19,282 pairs of ensemble forecasts and observations were available on 
291 distinct days and at 79 distinct observation locations. Additional data from the years 2006 and 
2007 were used to provide an appropriate rolling training period for all days in 2008, for the first 
phase estimation of the correlation model © and to establish the optimal length of the rolling train- 
ing period. Further information about the UWME, now using the WRF mesoscale model, as well as 



real time forecasts and observations, can be found online at |http : / /www . atmos . Washington 
| . edu/~ens/uwme . cgi} 



3.2 Training periods for Regional and Local EMOS 

As noted, we distinguish Regional and Local EMOS forecasts. The Regional EMOS method uses 
all available training data to estimate a single set of parameters that is used throughout the Pacific 
Northwest domain, and can be used directly on the model grid as well. The Local EMOS technique 
uses training data from the station at hand only to obtain a distinct set of parameters at each station. 
Thus, the method applies at observation stations only, and can not be used directly on the model 
grid. 

To determine a suitable value of the length n of the rolling training period for phases two and 
three of our estimation scheme, as described in Section [21 we considered experiments with wind 
vector forecasts based on the UWME in calendar year 2007. In these experiments, phase one of 
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Figure 6: Predictive performance of the Local EMOS method in calendar year 2007 as a function of the 
length of the rolling training period in terms of the relative mean energy score (left) and the relative bivariate 
mean absolute error (right). Each curve corresponds to an observation station in Washington state, with the 
dot indicating the lowest value of the performance measure, relative to the average value among the lengths 
considered. 



the estimation scheme used (here, in-sample) data from 2007, while phases two and three were 
done on a rolling training period of length n days. The predictive performance was evaluated 
using the mean energy score and the bivariate mean absolute error, as described in the Appendix. 
For the Regional EMOS method, these metrics differed by less than a half percent as the length 
n of the rolling training period varied betw een 25 and 60. Figure [6] summariz es the results for 
the Local EMOS method, where we follow iThorarinsdottir and Gneiting d2010l) and consider the 
observational locations in Washington state. Based on these results, our case study in calendar 2008 
uses a rolling training period of length n = 30 days for the Regional EMOS method, and of length 
n = 40 days for the Local EMOS technique. However, training periods of any length between 
n = 20 and n = 80 days work well, and the predictive performance of the EMOS methods is 
insensitive to this choice. 



3.3 Reference forecasts 

We compare the postprocessed Regional and Local EMOS forecasts to the raw UWME forecast, 
as well as to postprocessed reference forecasts, as described now. For all reference forecasts, we 
distinguish Regional and Local methods, using rolling training periods of the most recent available 
n = 30 and n = 40 days, respectively. 

A natural reference standard is the Independent EMOS technique, whi c h app lies the standard 
EMOS or heterogeneous Gaussian regression technique of iGneiting et alj (120051) to each of the 
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vector wind components u and v individually, and then combines them under the assumption of 
independence. This results in a postprocessed bivariate normal density forecast of the form (OQ) 
with the critical correlation parameter p uv constrained to be zero, but with the means and variances 
for u and v being essentially identical to those in the bivariate EMOS method, except for minor 
differences due to the slightly differing estimation schemes. In particular, the bivariate EMOS 
method and the Independent EMOS technique yield essentially identical deterministic, bivariate 
mean and median wind vector forecasts, even though the respective bivariate predictive densities 
may differ substantially. 



The Ensemble Copula Cou pling (ECC) method, originally hinted at byferemnes] ( 2007 ) and 
Krzysztofowicz and Tothl (|2008|) and recently investigated and developed by ISchefzikl (|2011l) . is a 



tool for restoring a raw ensemble's flow-dependent rank dependence structure in individually post- 
processed predictive distributions. Here, we describe and apply the method in the context of the 
Independent EMOS technique and a raw ensemble with m members. At any given location and pre- 
diction horizon, let u*, . . . , and v{, . . . denote samples from the individually postprocessed 
univariate EMOS distributions. As the postprocessed distributions for the wind vector components 
are univariate Gaussian, these are simply random numbers drawn from a univariate normal density. 
Let ui, . . . ,u m and v 1 , . . . ,v m denote the raw ensemble values for the respective wind components, 
and let r u and r v be permutations of the numbers 1 , . . . , m such that 

u Tu (i) < u Tu{2 ) < < u Tu{m) and v Tvil) < v Tv{2 ) < ■ ■ ■ < v Tv{m) , 
with any ties resolved at random. The ECC ensemble then consists of the m wind vectors 



Thus, the ECC ensemble inherits and honors the raw ensemble's rank dependence structure. For 
example, if the first raw ensemble member shows the second lowest u component and the third 
highest v component among the raw ensemble members, then the same property holds true for the 
EC C ensemb le. 

Schefzikl (2011) provides a detailed discussion of the ECC technique and exposes its ties to 



copulas dScholzel and Friederichs . 2008 ). The wind vector postprocessing technique recently pro- 
posed by IPinsonl (|201lh can be interpreted as a particularly attractive variant, where the values 



U 



u* are constructed as a translation and dilation of the raw ensemble values u±, 



u r 



and 



«*,...,t)^asa translation and dilation of v%, . . . , v m , respectively, rather than being drawn at ran- 
dom. Consequently, the postprocessed ensemble forecast retains both the raw ensemble's bivariate 
Spearman rank correlation and its (standard) Pearson product moment correlation structure. 

Finally, we consi der an Error Dressing ense mble as proposed by iGneiting et al.1 (|2008|) in the 
spirit of the work of iRoulston and Smith! (|2003|) . where we dress the UWME mean forecast with 
35 error vectors from the corresponding training set. 



3.4 Results over the Pacific Northwest 



We now give verification results aggregated over our test period, the calendar year 2008, and the 
Pacific Northwest domain of the UWME. For details about the verification techn iques method, 
which include the multivariate rank histogram, as proposed by IGneiting et al.l(|2008l) for calibration 
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Figure 7: Multivariate rank histograms for raw and postprocessed ensemble forecasts of surface wind vectors, 
aggregated over calendar year 2008 and the Pacific Northwest. 
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Table 2: Predictive performance of forecasts of surface wind vectors in terms of the mean energy score (ES) 
and the bivariate mean absolute error (bMAE), both in meters per second, and the reliability index A for the 
multivariate rank histogram, aggregated over calendar year 2008 and the Pacific Northwest. 



Method 


ES 


bMAE 


A 


Raw Ensemble 


LA 1 


J.Uf 




Regional Independent EMOS 


2.43 


2.79 


0.37 


Regional ECC 


2.33 


3.00 


0.02 


Regional Error Dressing 


2.19 


3.01 


0.01 


Regional EMOS 


2.01 


2.80 


0.03 


Local Independent EMOS 


2.28 


2.60 


0.21 


Local ECC 


2.16 


2.78 


0.07 


Local Error Dressing 


2.07 


2.87 


0.01 


Local EMOS 


1.87 


2.61 


0.03 



checks, and proper scoring rules (|Gneiting and Rafteryi . l2007f) . are given in the Appendix. All 
scores used are negatively oriented, that is, the smaller the better. 

Figure [7] shows multivariate rank histograms for the raw ensemble forecast and Regional and 
Local versions of the Independent EMOS, ECC, Error Dressing and (bivariate) EMOS techniques. 
The raw ensemble forecast shows a U-shaped rank histogram, thus indicating underdispersion. In 
contrast, the Independent EMOS forecasts are overdispersed, as they neglect to take dependencies 
between the wind vector components into account, which can be corrected for with the ECC tech- 
nique. In addition to the ECC method, the Error Dressing and (bivariate) EMOS techniques show 
uniform rank histograms, as expected from a calibrated forecast. 

Table[2]provides numerical summary measures of the predictive performance, with the reliabil- 
ity index A for the multivariate rank histogram confirming the visual impression in Figure [71 The 
energy score is a direct analogue of the continuous ranked probability score for univariate quanti- 
ties that provides an overall assessment of the quality of a probabilistic forecast, addressing both 
calibration and sharpness. The bivariate absolute error generalizes the absolute error and assesses 
deterministic forecast skill. 

Not surprisingly, the Local approaches outperform the Regional approaches, both in terms of 
the mean energy score and the bivariate mean absolute error. Also, the postprocessed Indepen- 
dent EMOS forecasts show lower scores than the raw forecast. The Independent EMOS and the 
(bivariate) EMOS forecasts have nearly identical bivariate mean absolute error, at a substantially 
lower value than for the other types of forecasts, including the ECC approach. This effect can be 
attributed to a discretization effect, in that the ECC technique turns the Independent EMOS density 
forecast into a discrete ensemble forecast. However, the ECC approach improves on the Indepen- 
dent EMOS forecast in terms of the energy score, as it honors the raw ensemble's flow-dependent 
bivariate dependence structure. The Error Dressing technique also shows good probabilistic fore- 
cast skill, as evidenced by a low energy score, even though it is unable to match the scores of the 
(bivariate) EMOS method, which performs the best. When compared to the raw ensemble, the 
postprocessed Local EMOS forecast reduces the mean energy score from 2.47 to 1.87 meters per 
second. 
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3.5 Results at Sea-Tac Airport 



Next we consider forecasts at the observation station at Sea-Tac Airport, with summary measures 
of the predictive performance in calendar year 2008 being provided in Table [3] The results mimic 
those for the Pacific Northwest. The Local EMOS forecast shows the highest probabilistic forecast 
skill, as quantified by the energy score, which decreases to 1 .94 meters per second, as compared to 
2.25 meters per second for the raw ensemble. 

For an illustration and explanation of how and to what extent the Local EMOS technique suc- 
ceeds in improving the raw ensemble forecast at Sea-Tac Airport, consider t he display in F igure [H 



which is a bivariate variant of the marginal calibration diagram proposed by iGneiting et al.l (|2007). 
Essentially, a marginal calibration diagram compares the observed climatology to the climatology 
incurred by the forecasts. In our display we plot, for each day for which data are available, the 
observed wind vector, perturbed very slightly in order to honor the undiscretized distribution and 
improve readability, a randomly chosen member of the UWME raw ensemble, and a wind vector 
sampled from the postprocessed, bivariate normal Local EMOS density forecast. Again we ob- 
serve that the raw ensemble is underdispersive and fails to predict any extreme wind vectors. The 
postprocessed Local EMOS forecast corrects for the underdispersion, thus leading to substantially 
improved probabilistic forecast skill, as evidenced by the energy score. However, the raw ensemble 
does not show any recognizable biases, and so the improvement in deterministic forecast skill is 
minor. 



3.6 Results for wind speed 

In addition to producing calibrated and sharp forecasts of wind vectors, the bivariate EMOS method 
can be used to predict wind speeds, and so can be compared to the technique of Thorarinsdottir and 
Gneiting (2010), which is custom tailored to this task. As wind speed is a nonnegative quantity, 
we employ truncated normal predictive distributions, using an estimation scheme that is based on 
optimum score estimation, both in Regional and Local versions, where we use rolling training 
periods comprising the most rece nt n = 30 and n = 40 days avail able, respectively. In what 



follows, we refer to the method of iThorarinsdottir and Gneitingl (|2010 ) as wind speed EMOS. 



To generate probabilistic forecasts of wind speed from the bivariate EMOS forecast, we sample 
one hundred wind vectors from the bivariate predictive distribution, and compute the Euclidean 
norm of each vector, thereby obtaining a discrete forecast ensemble of size m = 100 for wind 
speed. Table |4] compares the predictive performance of this approach to that of the wind speed 
EMOS technique. As noted, the wind speed observations are strongly discrete, with wind speeds 
below two knots recorded as zero, which applies to about 14% of the observations in the test period. 
From the perspective of wind vectors, these observations tend to fall into the center of the respec- 
tive predictive distribution, and so the effect of the discretization is weak. From the perspective of 
wind speeds, an observed value of zero is right at the boundary of the climatological range, and 
so the effect is nonnegligible. To account for it, we replace every observation of zero by a num- 
ber drawn uniformly and at random between zero and two knots, or between zero and 1.03 meters 
per seconds, when computing the corresponding continuous ranked probability score or absolute 
error. The bivariate EMOS technique, particularly in its Local versio n, nearly matches the predic 



tive p erformance of the specialized wind speed EMOS technique of IThorarinsdottir and Gneiting 



d2010b . 
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Table 3: Predictive performance of forecasts of surface wind vectors at Sea-Tac Airport in calendar year 
2008 in terms of the mean energy score (ES) and the bivariate mean absolute error (bMAE), both in meters 
per second, and the reliability index A for the multivariate rank histogram. 



Method 


ES 


bMAE 


A 


Raw Ensemble 


2.25 


2.77 


0.54 


Regional Independent EMOS 


2.41 


2.89 


0.29 


Regional ECC 


2.37 


3.06 


0.23 


Regional Error Dressing 


2.12 


2.92 


0.07 


Regional EMOS 


2.06 


2.90 


0.24 


Local Independent EMOS 


2.40 


2.73 


0.22 


Local ECC 


2.21 


2.87 


0.22 


Local Error Dressing 


1.98 


2.75 


0.02 


Local EMOS 


1.94 


2.74 


0.09 
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Figure 8: Marginal calibration diagram for forecasts of surface wind vectors at Sea-Tac Airport in calendar 
year 2008. For each day available, the plot shows the observed wind vector, a randomly chosen member of 
the raw ensemble, and a wind vector sampled from the postprocessed Local EMOS density forecast. 
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Table 4: Predictive performance of forecasts of wind speed in terms of the mean continuous ranked proba- 
bility score (CRPS) and the mean absolute error (MAE), in meters per second, averaged over calendar year 
2008 and the Pacific Northwest. 



Method 


CRPS 


MAE 


Raw Ensemble 


1.34 


1.68 


Regional Wind Speed EMOS 


1.11 


1.57 


Regional EMOS 


1.15 


1.63 


Local Wind Speed EMOS 


1.08 


1.50 


Local EMOS 


1.07 


1.51 



4 Discussion 

In this paper, we have proposed a bivariate EMOS approach to the statistical postprocessing of en- 
semble forecasts of wind vectors that results in bivariate normal density forecasts. In experiments 
with 48-hour ahead forecasts of surface wind vectors over the Pacific Northwest, based on the 
University of Washington Mesoscale Ensemble (UWME), the postprocessed EMOS density fore- 
cast proved to bias correct and calibrate the raw ensemble forecast, therefore resulting in strongly 
improved deterministic and probabilistic predictive skill. When compared to the raw ensemble 
forecast and aggregated over calendar year 2008 and the Pacific Northwest, the Local EMOS tech- 
nique reduced the bivariate mean absolute error by 13% and the bivariate mean energy score by 
24%. 

There are several directions into which our bivariate EMOS method could be develope d. In 



phases two and three of the estimation scheme, the exponential forgetting approach as in IPinson 



(|2011l) could be implemented, where the parameter estimates are updated in a computat i onally 



efficient, adaptive way. Furthermore, a geost atistical approach su ch as that of iKleiber et al.1 (1201 1|) 



in the context of Bayesian model averaging (R aftery et all |2005|) could be developed in order to 



spread the Local EMOS estimates, which currently are available at observation locations on ly, over 



the m odel grid. Alternatively, if an analysis is used to fit the EMOS model, as was done by IPinson 



(|2011|) . training data are available on the analysis grid, thereby allowing for a gridded Local EMOS 



approach. 

A closer look at the scatterplots in Figure [3] suggests that the conditional distribution of the 
observed wind vector, given the ensemble forecast, tends to be skewed. Thus, our bivariate EMOS 
approach, which currently is based on bivariate normal densities, could be extended to allow for 
skewed distributions, such as biv ariate skew-normal or ske w-t densities, similar to the wind vector 
time series methods proposed by lHering and Gentonl(|2010l) . However, any such approach would be 
considerably more complex, and as it is more difficult to estimate a more complex predictive model, 
it i s not clear whether or not such an extension would result in improved forecast performance. 

Sloughtend2009h and lSloughter et al.l (1201 1|) proposed a bivariate version of the Bayesian model 



averaging technique (BMA; Raftery et al., 2005) for postprocessing ensemble forecasts of wind 
vectors. Like our EMOS method, the BMA approach results in a bivariate forecast density. How- 
ever, the BMA forecast density is a finite mixture of bivariate, power-transformed normal densi- 
ties and thus can be multimodal, as opposed to the EMOS forecast density, which is necessarily 
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unimodal and elliptically symmetric. Thus, the basic setting of the EMOS approach is more par- 
simonious, thereby allowing for the key innovation in our work, namely the explicit modeling of 
the correlation between the u and v wind vector components, conditionally on the direction of the 
ensemble mean vector. 

In contrast to the BMA and EMOS methods, the ensem ble copula co upling technique (ECC; 
Schefzik 201 1) and the postprocessing approach proposed by lPinsonl (|201ll) generate discrete fore- 
cast ensembles that are constrained to have the same number of members, and the same bivariate 
Spearman rank correlation coefficients, as the raw ensemble. Therefore, these methods are par- 
ticularly well adapted to large ensembles, where there is no pronounced need for the transition 
from a discrete forecast ensemble to a density forecast, nor any need for statistical correctio n of the 



condit iona l correlation stru cture between the wind vector components. Accordingly, both iPinson 



d201lh and lSchefzikl (l201lh tailored their methods to the 50-member European Centre for Medium 
range Weather Forecasts (ECMWF) ensemble, while our work considered the much smaller eight- 
member UWME, where the wind vector ensemble forecasts for any given location and valid time 
may show unreliable, physically unrealistic empirical correlation coefficients. This effect is caused 
by the small ensemble size and corroborates the need for a parametric correlation model. 

In this paper, we considered probabilistic forecasts of wind vectors for a single location and 
valid time. While we modeled the bivariate correlation structure between the wind vector compo- 
nents, we considered the postprocessed bivariate density forecasts at each site individually, without 
modeling the dependence structure between locations. To address the latter, the methods devel- 



oped i n our paper could be combined with the spatial statist ical technique s introduced by lGel et al. 
(200 4|). A comm on advantage of the method proposed by IPinsonl (|2011[) and the ECC technique 
(ISchefzikL 1201 lb lies in their immediate extension to the calibration of spatio-temporal trajectories, 
in that the raw ensemble's bivariate rank dependence structure is inherited by the postprocessed 
ensemble, subject to the aforementioned caveats. 



Appendix: Verification methods 



In verifying probabil istic forecasts of a mu l tivariate weath er quant ity, we use techniques introd uced, 
studied and used by iGneiting et al.l (|2008|) . iGneitingl (|201 II) and lPinson and Hagedornl (|201 1|) . 

To assess the calibration of probabilistic forecasts of wind vectors, we use the multivariate rank 
histogram, which is a natural, direct generalization o f the verification rank histog r am or Talagrand 



diagram for a univariate quantity (|Andersonl 1 19961 : lHamill and Colucci. 1 19971 : iTalagrand et al. . 



19971) and can be interpreted analogously, in ways described by Ha milll (|200l|) . In particular, U- 
shaped multivariate rank histograms correspond to underdispersed ensembles, while inverse U- or 
hump-shaped histograms indicate overdispersed ensembles. For a calibrated ensemble, we expect 
a uniform rank histogram. For an ensemble with m members, the m ultivariate verificatio n rank is 
a possibly randomized number between 1 and m + 1, and we refer to IGneiting et al.1 (|2008|) for the 
technical details in its construction. To quantify the departure of the rank histogram from unifor- 
mity, we use the discrepancy or reliability index A proposed by lDelle Monache et al.1 (120061) . given 
by 

m+1 1 

E fi~ — TT > (6) 



A 
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where f\ is the observed relative frequency of verification rank i — 1, . . . , m + 1. The UWME raw 
ensemble and ensemble copula coupling (ECC) techniques result in a discrete forecast ensemble 
with m = 8 members, from which the construction of the multivariate rank histogram is straightfor- 
ward. The Error Dressing techniques provide discrete forecast ensembles with m = 35 members, 
and we bin the corresponding 36 ranks in the multivariate rank histogram into nine groups, com- 
prising ranks 1-4, . . . , 33-36, respectively, to facilitate the comparison. For the Independent EMOS 
and bivariate EMOS techniques, we draw a simple random sample of size m = 8 from the bivariate 
normal predictive distribution, and then compute the multivariate rank histogram. 

To assess the overall quality of p robabilistic forecasts ; considering both cal ibration and sharp- 
ness, we use proper scoring rules (IGneiting and Rafteryl . 120071: IWilksL 1201 1|) . For a univariate 
weather quantity, such as wind speed, the proper continuous ranked probability score is defined as 



crps(P, y) 



1 



[P{x) - l(x > y)f dx = E P \X -y\- -E P \X - X' 



(7) 



where P is the predictive distribution, here taking the form of a cumulative distribution function, 
X and X' are in dependent random variables with cumulative distribution function P, and y is the 
verifying value (IGneiting and Rafteryl . 120071) . The term I(x > y) denotes an indicator function, 
equal to 1 if x > y, and equal to otherwise, and E is the expectation operator. The absolute error 
is defined as 

ae(P,y) = \med P - y\, (8) 



where medp is a median of the probability distribution P (jGneitingL 1201 ll : IPinson and Hagedorn . 
l20Tlh . 

The energy score was introduced by Gneiting and Rafteryl ( 2007 ) and Gneiting et al. ( 2008 ) as a 
direct generalization of the continuous ranked probability score © in the evaluation of probabilistic 
forecasts of multivariate quantities. As we are interested in wind vectors, we restrict the discussion 
to bivariate weather quantities. The energy score then is defined as 



es(P,y) = E P \\X -y\\ - -E P \\X - X'\ 



(9) 



where || ■ || denotes the Euclidean norm in M 2 , P is the predictive distribution, X and X' are 
independent random vectors with distribution P, and y e M 2 is the verifying wind vector. For 
an ensemble forecast, the predictive distribution P cns has point mass ^ at the member forecasts 

Xi, 



-, and the energy score can be evaluated as 



^ m 

es(P cns ,y) 

m L — ' 



R - y\ 



1 

2m 2 



EE 

i=l j=l 



We use this formula to compute the energy score for the UWME raw ensemble, Ensemble Copula 
Coupling (ECC) and Error Dressing techniques. For the Independent EMOS and bivariate EMOS 



techniques, we draw a simple random sample xi, 



,x k e 



from the corresponding predictive 



density, and replace the exact energy score © by the computationally efficient approximation 



1 k 



1 fe-i 



2(Jfe-l) 



1/ 



X 



i+l 
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where we use a sample of size k = 10, 000. Similar approximations apply to the continuous ranked 
probability score 0. 

The natural generalization of the absolute error ([8]) is the bivariate absolute error 

bae(P, y) = ||bmedp — y\\, (10) 

where bmedp denotes the bivariate or spatial median of the probability distribution P, defined as 

bmedp = argmin xe]R 2 Ep\\x — X\\, 



where X is a random vector with distribution P (IVardi and ZhangL 120001 : iGneitingL 1201 1|) . For an 
elliptically symmetric distribution, such as a bivariate normal distribution, the bivariate median and 
the mean vector coincide. For other types of bivariate distributions, such as an ensemble forecast 



with point mass — at the member forecasts x%, 



the bivariate median generally is 



different from the corresponding mean vector, an d typically it need s to be determined numerically. 
For doing this we use the algorithm described by IVardi and Zhang fcOOOh and implemented in the 
R package ICSNP. 

In practice, forecasting methods are assessed by averaging scores over a test period, resulting in 
the mean continuous ranked probability score (CRPS), mean absolute error (MAE), mean energy 
score (ES) and bivariate mean absolute error (bMAE), respectively. All these quantities are nega- 
tively oriented, that is, the smaller the better, and we report their values in the unit of meters per 
second. 
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